Comments on tails in Schwarzschild spacetimes. 
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We performed a careful numerical analysis of the late tail 
behaviour of waves propagating in the Schwarzschild space- 
time. Specifically the scalar monopole, the electromagnetic 
dipole and the gravitational axial quadrupole waves have been 
investigated. The obtained results agree with a falloff l/t 2l+3 
for the general initial data and l/t 2l+4 for the initially static 
data. 



I. INTRODUCTION 

Waves propagating in a curved spacetime usually un- 
dergo backscatter; that leads to the emergence of two 
classes of effects, the so-called quasinormal modes (dis- 
covered by Vishveshvara [1]) and the tails. Price, who 
investigated the temporal behaviour of tails in 1972, has 
shown some kind of universality. The late tail behaviour 
of scalar, electromagnetic and gravitational waves hap- 
pened to depend on the order of the multipole expansion 
[2] . The exposition of [2] is somewhat confusing and that 
forces us to say more than would normally be needed for 
the purpose of a paper reporting only numerical results. 

Price states If a static l-pole field is present outside the 
star, prior to the onset of collapse, the field will fall off as 
t~ ( 2l + 2 ) . If there is no field initially outside the star, but 
an l-pole perturbation develops during the collapse pro- 
cess, it will fall off as £-( 2i + 3 ) a t large time [2]. That re- 
quires a comment. Since the evolution equation (2) (see 
below) is linear and of the form — d 2 ^ + L^> = 0, where L 
is a linear time-independent operator, its solutions can be 
superposed. That means that if there is a static solution, 
then it cannot influence any evolving perturbation. Irre- 
spective of whether there are static solutions or not, the 
perturbation will decay in the same way. Therefore the 
quoted statement cannot be taken at the face value. The 
natural understanding would be the following: A static 
field means that there is initially a field = 0) ^ 0), 
but dt^ — (hence data are momentarily static). The 
initial absence of field in turn means that 'J = but 
9f\& 7^ 0. This agrees with the standard way of posing 
the initial value problem in mathematical literature, and 
with this in mind we inspect the later papers on this sub- 
ject. The notion of initially static l-pole has been clearly 
specified in two publications that appeared in late seven- 
ties; it becomes an initially stationary multipole of order I 
(that falls off like l/t 2l+2 ) [3]. Therein the term initially 
stationary means momentarily static as defined above, 
but with \& satisfying special asymptotic condition. On 



the other hand a detailed discussion in [4] shows that for 
initially nonstatic, dt^ ^ analytic initial data satisfy- 
ing certain boundary conditions (at the event horizon of a 
black hole and at spatial infinity) the decay of a solution 
agrees with that of the Price's initially static l-pole field. 
This might suggest in turn that one can understand the 
" static l-pole field" as a specially prepared general initial 
wave profile. Corresponding initial data (both \& and dt^ 
are nonzero) have a noncompact support (the spatial sup- 
port has infinite extension, as measured in terms of the 
tortoise variable r* - see below). The work of Gundlach 
et al. [5] numerically confirmed the analytically derived 
l/t 2l+2 decay of late tails. Let us remark an internal in- 
consistency - the contention must be now that the fall off 
actually depends on the profile of initial data, in contrast 
to the primary expectation that the "tail" is a universal 
phenomenon due to the backscatter of the wave signal off 
the curvature of the geometry. 

Leaving aside the area of initially static fields, let us 
notice that in the case of data of compact support the 
predictions of [2] are clear - solutions should have a falloff 
£-(22+3)^ xhat estimate applies to solutions generated 
by generic initial data ^ 0, / 0) [4]. From a 
naive inspection of the formal solution \1/ = dtG * ^(t = 
0) + G * dt^it — 0) (where G is the Green function), one 
expects that in the case of momentarily static initial data 
the fall off is faster (by 1/t) than that corresponding to 
general initial data. It is interesting that no numerical 
investigation has been performed in the case of moment 
of time symmetry initial data of compact support. It 
should be noticed also that the status of the analytic in- 
vestigation is still not quite satisfactory, since it consists 
either of heuristic analysis [2] or rather formal investiga- 
tion of the Green function formulae ( [4], [6]) [8]. It is 
only recently that the (generic case) falloff 1 /t 3 has been 
rigorously proven by Dafermos and Rodnianski for the 
monopole mode of the real scalar field [9] . The same con- 
clusion is being derived, in a different (spectral) method, 
by Machedon and Stalker [10]. 

The main focus of our numerical work will be on find- 
ing the falloff of tails generated by initial data of compact 
support in either of the aforementioned two fundamen- 
tal cases. We study the evolution of the scalar monopole 
(1 = 0), the electromagnetic dipole (I = 1) and the grav- 
itational axial quadrupole (/ = 2). The obtained results 
reveal a falloff l/t 2l+3 for the data d t ^i ^ 0, = 0, 
in agreement with [2], and l/t 2l+i for the initially static 
initial data [11]. This latter result on the behaviour of 
late tails generated by initial data d t ^>i = 0, ^ is 
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new in the numerical literature. Thus the generic initial 
data have the falloff l/t 2l+3 , which agrees with Leaver 
[4], Ching et al. [6] and Poisson [7]. 

Our results support the view that the 1-th moment 
of any of the waves (scalar, electromagnetic or gravita- 
tional) will decay like l/t 2l+3 for general initial data (in 
accordance with [6], [4] and [7]) and l/t 2l+4 for the 1-th 
moments evolving from initially static data (in agreement 
with the formal analysis of [4] and [7]). These results 
should be of practical significance for numerical relativ- 
ity. The determination of the late tail behaviour is a 
nontrivial but at the same time feasible numerical task 
that can serve as a useful test for the accuracy and sta- 
bility of numerical codes. 



II. DEFINITIONS AND EQUATIONS 

The spacetime geometry is defined by the line element 

,/G = -1— + /?-(/()-. (1) 



Jn2 

- m dt 2 + — + R 2 dn 2 , 
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where t is a time coordinate, R is the radial areal coordi- 
nate, 7]r = 1 - and dfl 2 = d9 2 + sin 2 6d<j) 2 is the line 
element on the unit sphere, < <fi < 2ir and < 9 < it. 
Throughout this paper the Newtonian constant G and 
the velocity of light c are put equal to 1. 

The propagation of the scalar, the dipole electromag- 
netic and the axial (quadrupole) gravitational waves is 
given by 



{-d 2 + d%)^> = V*. 



(2) 



Here r*(R) = R + 2mln^2^; — l) is the tortoise coordi- 
nate while the potential term reads: for the I = mode 
of the scalar field 



V(R) = 2m^, 
for the I = 1 (dipole) electromagnetic mode 

V(R) = 2^ 
and for the quadrupole axial mode of GW 



(3) 



(4) 



(5) 



III. NUMERICAL SCHEMES 

We will choose following classes of initial data: 

i) *(R,t = 0) = sin 4 (^r^M) and ^| t=0 = 

-<9 r .\E' 4= o for r* G (r*(a) = r*(3m),r*(a) + 4m), and 
*(r, 0) = = for r* < r*(a) or r* > r*(a) + 4m. 



The expected asymptotic tail falloff should be the same 
as for generic initial data. 

ii) Symmetric initial data with \P(i?, t = 0) = 

sin 4 (7T r ^ (a) ) for r* G (r*(a) = r* (3m), r* (a) + 4m) 
and *(r, 0) = for r* < r*(a) or r* > r*(a) + 4m, and 
dt^t=a = everywhere. It is shown below that the cor- 
responding tail decays faster than that of i). 

iii) d t V = sin 4 ( V*- 4 ^ (a) ) for r* < r*{a) = r*(3m) 

or r* > r*(a) + 4m and d t ^(r,0) = for r* < r*(a) or 
r* > r*(a) + 4m, while \&|t=o = everywhere; the tail 
behaviour happens to be like in i). 

One finds convenient to split the wave equation (2), 
into the pair of first order differential equations 



(d t + d r »)* = $ 



(6) 



(d t -d r *)$ = -V(r*)<H. (7) 
We use the following mesh scheme (see Fig. 1) 



*ab 



dt 



a1 



*b1 



*a2 



*b2 



dr = 2 dt 

FIG. 1. This figure shows the mesh used in the numerical 
computations. In the one-step scheme the functions $ a f>, ^G;, 
are computed from Eqs (9 - 11) using values at points al and 
61. The second algorithm needs values at four points al, a2, bl 
and 62 in order to obtain from Eqs (14). 

and two differencing schemes: the one-step implicit al- 
ghoritm and the multistep Adams-like algorithm. In the 
first approach we approximate the equations (7) by 



dt 



(8) 



dt 



®ab - $61 = - (Vab^ab + V bl * 61 ) — . (9) 

In the next step these equations are solved for ^ a b, & a b, 

*a6 = *a&(*al, * 61, $al, *6l) (10) 



$a6 = $o&(*ol, $al, $&l)- 



(11) 



The second alghoritm is based on the Adams-Moulton 
corrector formula for ordinary differential equations [12] 

Vn+i = Vn + y^(%„+i + 8y„ - y' n _{) + 0{dt A ) (12) 

Applying this formula to each of the equations (7) we get 
approximately 
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*a6 = *al + ^(5$ a6 + 8$ a i-$ a2 ) (13) 

and 

$ah = $61 - ^(5V„6*a6 + 8V 6 1*61 - H 2 *w) (14) 

Again this linear set of equations can be solved for 

*ab = *ofc(*ol, *61, *o2, *62, $al, $61, *a2, $62) 

$afc = *o&(*al, *61, *o2, *62, $al > $61 , $a2, $& 2 )- 

(15) 

Let us stress that these two alghoritms give almost the 
same results, hinting at the numerical stability of our 
methods. 



IV. NUMERICAL RESULTS 

We expect the field ^i(r*,t) to behave (for a fixed r*) 
like 



The numerical time ranged from 5000 m (for the grav- 
itational and electromagnetic waves) up to 40000 m (for 
the scalar waves) . It turned out that the precise shape of 
this function, without numerical noise, can be calculated 
using numbers with 64 significant digits (gravitational 
and electromagnetic cases) and the ordinary double pre- 
cision numbers (scalar case). In order to achieve this 
aim we have used the freely distributable "qd" (quad 
precision) and "arprec" (arbitrary precision) numerical 
libraries [13]. 




lim tf(r*,t) = Ct~ 



(16) 



1000 2000 3000 4000 5000 

FIG. 3. The solid (general initial data) and broken line 
(data with = 0), the behaviour of f(t) for the selectromag- 
netic waves. 



Our goal is to get a value of a with a reasonable accuracy. 
Therefore we calculate the function 



/(*) 



dlog(tt t (r*,t)) 

dlog(t) 



-t 



rflog(tt t (rV)) 
dt 



(17) 



which asymptotically should be equal to a. 

In the first instance we analysed the case with waves, 
which belongs to the class i) of Sec. 3, and for the initial 
data of class iii) of the preceding section. The results 
are presented in Figs. 2-4 which show the temporal 
behaviour of the funtion f(t) as seen at r* = r*(3m)+4ra. 
In all figures (2 - 7) below the time is put on the x-axis, 
which is scaled in units of m. 
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FIG. 2. The solid (general initial data) and broken line 
(data with = 0), the behaviour of f(t) for the scalar waves. 
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The solid (general initial data) and broken line 
= 0), the behaviour of f(t) for the gravitational 



The high accuracy calculations are time consuming 
and therefore the achieved calculational time for the 
scalar case exceeds by a factor of 8 the evolution time 
for the remaining cases. The numerically achieved values 
of the exponent , for the general initial data i), range 
from 2.999 for the scalar field, through 4.99 for the elec- 
tromagnetic field up to 6.99 for the gravitational waves 
(figures 2-4, solid line). Both sets of exponents is in a 
good agreement with 3, 5 and 7, respectively, obtained 
by Ching et al. [6]. Similar results (broken line Fig. 3) 
have been obtained for the initial data of class iii) of the 
preceding section, in agreement with [2]. 

The momentarily static initial data (class ii) from Sec. 
3) produced profiles f(t) shown in Figs. 5-7. In this 
case the numerical exponents have been determined at 
r*(b) = r*(3m) + 4m and also at another obervation 
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point, r* = 100m + r*(b). In the first case we obtained 4, 
6 and 8 (up to the fourth digit number) while in the latter 
case we arrived at 3.99, 5.88 and 7.84, for the scalar, elec- 
tromagnetic and gravitational waves, respectively. These 
exponents seem to be stabilize at 4, 6 and 8 for the scalar, 
electromagnetic and gravitational fields, respectively; no- 
tice that the results detected at r* = r* (3m) + Am con- 
verge quicker to the prospective limiting value than those 
taken at the point located farther. The asymptotic tail 
regime is evidently achieved quicker at a region closer to 
the horizon, which is a new feature, unknown in the exist- 
ing literature. But in both detection points, the asymp- 
totic values of a are similar and very close in the case 
of the scalar waves, see Fig. 5. The calculation of these 
exponents constitutes the main result of this paper. 
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FIG. 5. Time-symmetric initial data. The solid and bro- 
ken lines show the behaviour of f(t) for scalar waves, as 
seen at the observation points r*(b) = r*(3m) + 4m and 
r* = r*(b) + 100m, respectively. 
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FIG. 6. Time-symmetric initial data. The solid and broken 
lines show the behaviour of f(t) for electromagnetic waves, 
as seen at the observation points r*(b) = r*(3m) + Am and 
r* = r*(b) + 100m, respectively. 
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FIG. 7. Time-symmetric initial data. The solid and bro- 
ken lines show the behaviour of fit) for gravitational waves, 
as seen at the observation points r*(b) = r*(3m) + 4m and 
r* = r*(b) + 100m, respectively. 
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